function a = sim_all(param,o)
global theta e12 eg ey;

gamma1=param(1);  gamma2=param(2);
betay=param(3);   cy=param(4);      sigtheta=param(5);    
c1=param(6);      c1d=param(7);     phi1=param(8);    phi1d=param(9);    beta1=param(10);  beta1d=param(11);
c2=param(12);     c2d=param(13);    phi2=param(14);   phi2d=param(15);   beta2=param(16);  beta2d=param(17);
cg=param(18);     cr=param(19);     cm=param(20);
phig=param(21);   phir=param(22);   phim=param(23);
sigg=param(24);   sigm=param(25);   sigr=param(26);
if o ==1 % other race
 %   gamma1 = gamma1+gamma12; gamma2 = gamma2+gamma22;
    c1 = c1+c1d; c2=c2+c2d;
    phi1 = phi1+phi1d; phi2 = phi2+phi2d;
    beta1=beta1+beta1d; beta2=beta2+beta2d;
end;

g=cg+phig*theta+eg;
Tstar = [c1+phi1*theta+beta1*g,c2+phi2*theta+beta2*g]+e12;
T1star = Tstar(:,1); T2star = Tstar(:,2);
T1=(T1star>=0)+0; T2=(T2star>=0)+0;



b1 = T1-normcdf(cy+theta+betay*g);
b2 = T2-normcdf(cy+theta+betay*g);

Ystar = cy+theta+b1*gamma1+b2*gamma2+betay*g+ey;
Y=(Ystar>=0)+0;
PrY = normcdf(cy+theta+b1*gamma1+b2*gamma2+betay*g);

a=[Y, T1, T2, b1, b2,T1star,T2star,Ystar,PrY];



